### texreg method for polr objects
extract.pooled.polr <- function(mod.list, include.thresholds = TRUE, include.aic = TRUE,
                         include.bic = TRUE, include.loglik = TRUE,
                         include.deviance = FALSE,
                         include.nobs = TRUE, include.lrtest = TRUE,
                         include.lrtest.divided = TRUE, ...) {

  lrtest.chi <- median(do.call("rbind", lapply(mod.list, function(x){round(lmtest::lrtest(x)$C[2], 3)})))
  lrtest.pval <- median(do.call("rbind", lapply(mod.list, function(x){round(lmtest::lrtest(x)$P[2], 3)})))

  lrtest.divided.chi <- median(do.call("rbind", lapply(mod.list, function(x){round(lrtest(x, . ~ . - dividedpressenate)$C[2], 3)})))
  lrtest.divided.pval <- median(do.call("rbind", lapply(mod.list, function(x){round(lrtest(x, . ~ . - dividedpressenate)$P[2], 3)})))
  
  coef.frame <- do.call("rbind", lapply(mod.list, coef))
  zeta.frame <- do.call("rbind", lapply(mod.list, function(x){x$zeta}))
  
  coefs <- apply(coef.frame,2,mean)
  ses <- apply(coef.frame,2,sd)
  zetas <- apply(zeta.frame, 2, mean)
  zeta.ses <- apply(zeta.frame, 2, sd)
  
  coef.names <- colnames(coef.frame)
  zeta.names <- colnames(zeta.frame)
  
  threshold.names <- zeta.names
  threshold.coef <- zetas
  threshold.se <- zeta.ses
  threshold.zval <- threshold.coef / threshold.se
  threshold.pval <- 2 * pnorm(abs(threshold.zval), lower.tail = FALSE)
  
  beta.names <- coef.names
  beta.coef <- coefs
  beta.se <- ses
  beta.zval <- beta.coef / beta.se
  beta.pval <- 2 * pnorm(abs(beta.zval), lower.tail = FALSE)
  
  if (include.thresholds == TRUE) {
    names <- c(beta.names, threshold.names)
    coef <- c(beta.coef, threshold.coef)
    se <- c(beta.se, threshold.se)
    pval <- c(beta.pval, threshold.pval)
  } else {
    names <- beta.names
    coef <- beta.coef
    se <- beta.se
    pval <- beta.pval
  }
  
  n <- mean(do.call("rbind", lapply(mod.list, nobs)))
  lik <- median(do.call("rbind", lapply(mod.list, logLik)))
  aic <- median(do.call("rbind", lapply(mod.list, AIC)))
  bic <- median(do.call("rbind", lapply(mod.list, BIC)))
  dev <- median(do.call("rbind", lapply(mod.list, deviance)))
  lrt <- lrtest.chi
  lrp <- lrtest.pval
  lrt.d <- lrtest.divided.chi
  lrp.d <- lrtest.divided.pval
  
  gof <- numeric()
  gof.names <- character()
  gof.decimal <- logical()
  if (include.aic == TRUE) {
    gof <- c(gof, aic)
    gof.names <- c(gof.names, "AIC")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.bic == TRUE) {
    gof <- c(gof, bic)
    gof.names <- c(gof.names, "BIC")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.loglik == TRUE) {
    gof <- c(gof, lik)
    gof.names <- c(gof.names, "Log Likelihood")
    gof.decimal <- c(gof.decimal, TRUE)
  }

    if (include.deviance == TRUE) {
    gof <- c(gof, dev)
    gof.names <- c(gof.names, "Deviance")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  
  if (include.lrtest == TRUE) {
    gof <- c(gof, lrt, lrp)
    gof.names <- c(gof.names, "Likelihood Ratio Test", "%Likelihood Ratio p-value")
    gof.decimal <- c(gof.decimal, TRUE, TRUE)
  }

  if (include.lrtest.divided == TRUE) {
    gof <- c(gof, lrt.d, lrp.d)
    gof.names <- c(gof.names, "Likelihood Ratio Test of Significance of Divided Government", "%Divided Likelihood Ratio p-value")
    gof.decimal <- c(gof.decimal, TRUE, TRUE)
  }

  if (include.nobs == TRUE) {
    gof <- c(gof, n)
    gof.names <- c(gof.names, "Number of Observations")
    gof.decimal <- c(gof.decimal, FALSE)
  }
  
  tr <- createTexreg(
    coef.names = names,
    coef = coef,
    se = se,
    pvalues = pval,
    gof.names = gof.names,
    gof = gof,
    gof.decimal = gof.decimal
  )
  return(tr)
}
